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Abstract. In this article the notion of metabolic turnover is revisited in the light 
of recent results of out-of-equilibrium thermodynamics. By means of Monte Carlo 
methods we perform an exact uniform sampling of the steady state fluxes in a genome 
scale metabolic network of E Coli from which we infer the metabolites turnover times. 
However the latter are inferred from net fluxes, and we argue that this approximation 
is not valid for enzymes working nearby thermodynamic equilibrium. We recalculate 
turnover times from total fluxes by performing an energy balance analysis of the 
network and recurring to the fluctuation theorem. We find in many cases values one 
of order of magnitude lower, implying a faster picture of intermediate metabolism. 
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The recent wealth of data coming from genome sequencing in biology|Xj is eager 
for unifying schemes, interpretations and insights that could come from physics. 
In particular metabolism, the ubiquitous and highly conserved enzymatic network 
devoted to free energy transduction in every cell, has been the subject of structural 
reconstructions at the scale of the whole genome [2]. Metabolism has deep physical 
roots and thus a long standing tradition of physical modelization effortsp]. A current 
challenge faced by physics thus concerns the extension of such efforts in large scale 
models. On one hand we lack detailed information on many parameters, on the other, 
even simple models with minimal assumptions lead to difficult computational issues. 
Much attention has been payed to the structural properties of metabolic networks pf], 
but on the other hand metabolism is inherently dynamical and a fundamental inherently 
physical question regards the assessment of its typical timescales. In this article we will 
analyze metabolites turnover times at the scale of the whole genome in a metabolic 
network model of E Coli. We will point out that such analysis requires to integrate 
thermodynamic information and thus the evaluation of an energy balance analysis of 
the network pi E]. 

If we model a metabolic system in terms of the dynamics of the concentration levels, 
assume well-mixing (no space) and neglect noise (continuum limit), we still have a 
very large non-linear dynamical system whose parameters can be unknown. For a 
chemical reaction network in which M metabolites participate in N reactions with the 
stoichiometry encoded in a matrix S = {S^r}, the concentrations c M change in time 
according to mass-balance equations 


c = S ■ v (1) 

where Vi is the flux of the reaction i that is in turn a function of the concentration levels 
Vi( c). Upon considering a steady state (homeostasis), i.e. a flux configuration satisfying 

S ■ v = 0 (2) 


we could in principle determine rigorously the timescales by performing a linear stability 
analysis of these steady states, i.e. upon linearizing the laws Vj(c) and finding the 
spectrum of the resulting matrix. Such calculation requires knowledge of the elasticity 
coefficients (7jj and in turn of the reaction laws with their parameters, that is not 
the case in large scale models. A widely employed approximation, if at least fluxes and 
concentrations are known, is to consider the metabolites turnover times[8],[9] t, i.e. the 
ratio between the concentration of a given metabolite c and the flux of production P 
(or equivalently destruction D } that is the same in the steady state), i.e. 


c = P — D = 0, 


c 

T= P 


(3) 


This turnover time is usually directly and intuitively interpreted as the average time 
it takes to fully replenish a given metabolic pool. However, we point out that, 
especially nearby thermodynamic equilibrium, net fluxes result from the difference 
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between forward and backward contributions v = u + — u~. This would imply that 
production and destruction fluxes split as well and the resulting turnover time could be 
lower [TO]: 


c = (P + + D~) - ( D + + P~) = 0 


(4) 


It should be understood that these contributions directly affect even the calculation 
of relaxation times from a more rigorous linear stability analysis. Now, if we have 
information about the net flux v and the free energy AG, we can estimate the backward 
and forward contribution v + ,v~ from a simplified form jllj of one of the main result in 
out-of-equilibrium thermodynamics, the fluctuation theorem[1211TBI 14]: 


— = e- AC/RT (5) 

v~ 

It has been shown that this simplified form is valid for the mass action law and the 
reversible Michellis-Menten kinetics, and it has been proposed to be valid in general [TT]. 
For instance, consider a metabolite in the red cell, Glucose-6-phosphate, produced 
in glycolysis na by an irreversible reaction, Hexokinase (AGi ~ —29 KJ/mol), and 
consumed by a reversible one, phosphoglucoisomerase (AG 2 — —2.9 KJ/mol). We can 
calculate at RT = 2.5KJ/mol that the turnover time To from net fluxes overestimates 
the one r from total fluxes by a factor 


Tq-t 

T 


AG, 

e rt 


~ 45% 


- 1 


( 6 ) 


In the following we will show the results of an uniform sampling of the possible steady 
state fluxes in a genome scale metabolic network model for the bacterium E Coli. Then, 
the results of an energy balance analysis of the network will be presented in order to 
estimate total fluxes from the simplified form of the aforementioned fluctuation theorem. 
We will calculate metabolites turnover times from net fluxes, correct them from estimate 
of total fluxes and show that in the latter case they can be much lower. We will finally 
draw out some conclusions, for instance regarding how this new faster picture of the 
intermediate metabolism affects the well-mixing hypothesis. 


2. Results and discussion 

We consider the steady state fluxes of the metabolic network model of E Coli iJR904[T6] 
in a glucose-limited minimal medium (see materials and methods sec.). 

In constraint based modeling, apart from mass balance constraints, fluxes are 
bounded in certain ranges v r G [n™ m ,n™ ax ] that take into account thermodynamical 
irreversibility, kinetic limits and physiological constraints. The set of constraints 


S ■ v = 0, 

V r G [u r min ,U r maX ] 


( 7 ) 





Metabolic turnover from EBA 


4 


defines a convex closed set in the space of reaction fluxes: a polytope from which feasible 
steady states can be efficiently inferred with Monte Carlo methods[17] (see materials 
and methods sec.) 

Once we have the flux distributions, we can single out for each metabolite the net 
production flux, that is the sum of positive definite terms (and that it is equal to the 
net consumption flux under our steady state assumption), if we have information about 
the concentration levels we can thus calculate the turnover times, i.e. 


— Pfx — 0 

Pfi = y ^OiSinV^SjuVi 

i 

Dfj, = - 


T = PE 
V p 


( 8 ) 


We show in fig 1 (top) in ordered fashion the mean of such production fluxes 
in log scale (for the sake of avoiding confusion we indicate the name of only some 
metabolites on the x axis). They span 8 orders of magnitude ranging from 10 1-2 to 10 - ' 
mmol/gDWh, thus giving an highly heterogeneous picture of intermediate metabolism. 
This is consistent with the time-hierarchy hypothesis, by which typical metabolic 
timescales should be highly heterogeneous in order to suppress instabilities [HI [19, [20] . 
As we mentioned in the introduction we can go a step beyond and calculate forward 
and backward contribution of reaction fluxes by using the fluctuation theorem from the 
knowledge of the reactions’ free energies. These can be estimated by performing an 
energy balance analysis of the network. Reactions free energies should be consistent 
with the reaction directions: sign(uj)AGj < 0 and if we decompose them in terms of 
metabolites chemical potentials A Gi = SyAA we have for given reaction directions 
the feasible space of chemical potentials: 


= sign (vi)S itl (9) 

'y ^ — o 

from which free energies can be inferred with relaxational methods starting from a 
reliable experimental prior. 

Once a free energy vector has been retrieved, the backward and forward 
contributions to fluxes can be calculated from the fluctuation theorem: 


Vi 

1 _ e A Gi/RT 
Vi 

e ~AGi/RT _ i 

( 10 ) 


The total fluxes are reported in fig 1 (bottom) in log scale comparing them with net 
values. 
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Figure 1. Top: Net metabolites production fluxes from uniform sampling of steady 
states of the metabolic network model of E Coli iJR904 in a minimal medium. Bottom: 
Total fluxes from energy balance analysis, compared with net fluxes. 


We can see that for many metabolites the corrected production flux can be one 
order of magnitude higher when it is processed by enzymes working nearby thermal 
equilibrium, that would imply faster relaxation times of the system, for example in its 
response to perturbations. In figure 2 (top) we plot for instance the distribution of ATP 
production fluxes in log scale: we pass from a distribution centered around the mean 
with value 56mmol/gDWh to an heterogeneous distribution with a fat tail and peaked 
around 350nnnol/gDWh. This correction lead to different qualitative estimate of the 
relative turnover time: in fig 2 (bottom) the distribution of the turnover time of ATP 
inferred from net and total fluxes assuming a concentration of 9.6mM|I?Tj is reported: we 
pass from a mean of 2.0s (consistent with previous estimates reported in databases[22]) 
to 0.4s upon using the correct value from total fluxes. 

In fig 3 we plot in log scale turnover times from net and total fluxes for several 
metabolites for which we have used the measures of concentration of|2T] We present in 
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Figure 2. Top: net and total ATP production flux distributions in the metabolic 
network model of E Coli iJR904 in a glucose-limited minimal medium. Bottom: 
inferred ATP turnover time distributions. 


Metabolite 

Tunover time from net flux (s) 

Turnover time from total flux (s) 

Glutammate 

90 ±20 

16 ±8 

3-Phosphoglycerate 

2.0 ±0.2 

0.12 ±0.6 

ATP 

2.0 ±0.1 

0.4 ±0.2 

ADP 

0.120 ±0.005 

0.02 ±0.01 

AMP 

0.5 ±0.1 

0.11 ±0.06 

NAD 

1.1 ±0.1 

0.3 ±0.1 

NADH 

3.5 ± 0.2 • 10 -2 

1.0 ±0.4 • 10' 2 

NADP 

i .6 ± 0 . 2 - nr 3 

2 ± 1 • 10 -4 

NADPH 

9 ±1•10“ 2 

2 ± 1 • 10“ 2 


Table 1. Turnover time estimates from net and total fluxes 
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Figure 3. Turnover time of several metabolites inferred from net and total fluxes. 

table 1 for comparison the turnover times estimate from net and total fluxes of some 
key metabolites. 

3. Materials and methods 

3.1. Data 

The network employed is E Coli iJR904pj]. This network consists of N = 1075 reactions 
among M = 761 metabolites. Upon considering a glucose limited minimal medium we 
are left, after removing blocked reactions and leaf metabolites with N = 667 reactions 
among M = 450 metabolites. We set a maximal glucose uptake of u g = —6mmol/gDWh 
and fix a minimal ATP maintenance of vatp = 7.6mmol/gDWh. The resulting polytope 
has dimension D = 233. In regard to the data for the energy balance analysis we recur 
to the chemical potential estimates coming from a genome scale application of the group 
contribution methods from [231124] , In order to calculate turnover times we have recurred 
to the measure of absolute concentrations reported in|21j. 

3.2. Computational methods 

3.2.1. Uniform sampling of steady states An uniform sampling of a convex polytope in 
high dimension is usually performed with Markov chain Monte Carlo methods, since an 
exact enumeration of the vertex would be infeasible due to their exponential number and 
static rejection methods[25j would suffer as well from high-dimensionality issues. We 
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have recurred to well known hit and run markov chain[26] that has been investigated 
in much detail in mathematical literature, showing nice convergence properties: it is 
guaranteed by detail balance to converge to the uniform distribution with a mixing 
time that scales as [27] 

t = o(d\ ( 11 ) 

where r and R are the radius of respectively the biggest (smallest) inscribed (inscribing) 
sphere. The factor R/r (’sandwitching ratio’) can lead to ill-conditioning issues, since 
the timescales of metabolic fluxes are typically very heterogeneous. Such factor can be 
reduced to a polynomial of the space dimension with a polynomial time algorithm that 
finds a rounding ellipsoid[25]. We have rounded the polytope with an ellipsoid founded 
with the flux variability analysis[29], that has been shown to reduce the sandwitching 
ratio to values that allow an efficient sampling|30j 

3.2.2. Energy balance analysis Finding a free energy vector consistent with a given flux 
configuration amounts at solving a system of linear inequalities for which relaxational 
algorithms can be employed .'IT. 32] starting from the experimental prior that we have 
reported in the subsec. Data. However we point out that duality arguments[33. [34] 
lead to the infeasiblity of such system if closed reactions loops are present. Such closed 
reaction loops has been exhaustively enumerated in this network and were corrected in 
a minimal way [35]. 

4. Conclusions 

In this work we have revisited the notion of metabolic turnover in the light of recent 
results in out-of-equilibrium thermodynamics. The net flux of a reaction working nearby 
thermodynamic equilibrium result from a contribution of backward and forward fluxes, 
whose value can be inferred from the fluctuation theorem upon knowledge of the free 
energy. The higher resulting total fluxes lead to effectively faster relaxation times in 
metabolic systems, a notion that can be captured at an approximated level by the 
computation of metabolites turnover times. We have performed an uniform sampling of 
the steady states of an E Coli genome scale constraint based metabolic network model, 
we have peformed an energy balance analysis of the network and we have estimated 
total production fluxes from the fluctuation theorem. We have shown that metabolites 
turnover times estimated in this way can be as far as one order of magnitude lower 
than the ones inferred from net fluxes. Such reduction of turnover times could lead in 
principle to values that are below typical metabolites diffusion times and thus it could 
affect the well-mixing hypothesis, that is at the core of our approach and constraint- 
based modeling in general. An approximate estimate[22] in which we consider a diffusion 
constant of D ~ 200 fim 2 / s, and that the diameter of E coli is d ~ 1 urn , lead to an order 
of magnitude estimate td — ^ — 10~ 3 s. The turnover time of all the metabolites we 
have calculated is above this threshold apart from NADP and adenosine. Such analysis 
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could be performed in principle for any metabolic network upon knowledge of the free 
energies landscape. On the other hand a more rigorous estimate of the true relaxation 
times of a metabolic system would requires genome-scale insights on enzymatic kinetic 
laws, including allosteric regulations, an aspect that would require further investigations. 
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